Method for Computing Short-Range Forces between Solid-Liquid Interfaces driving 

Grain Boundary Premelting 



00 

o 
o 

(N 
Oh' 

CO 



J. J. Hoyt 1 , David Olmsted 2 , Saryu Jindal 3 , Mark Asta 3 , Alain Karma 41 

Department of Materials Science and Engineering, McMaster University, 
Hamilton, ON, Canada 
2 'Sandia National Laboratories, Albuquerque, NM 
3 Department of Chemical Engineering and Materials Science, University of California, Davis, CA 
4 Department of Physics and Center for Interdisciplinary Research on Complex Systems, Northeastern University, Boston, MA 

We present a molecular dynamics based method for computing accurately short-range structural 
forces resulting from the overlap of spatially diffuse solid-liquid interfaces at wetted grain boundaries 
close to the melting point. The method is based on monitoring the fluctuations of the liquid layer 
width at different temperatures to extract the excess interfacial free-energy as a function of this 
width. The method is illustrated for a high energy E9 twist boundary in pure Ni. The short-range 
repulsion driving premelting is found to be dominant in comparison to long-range dispersion and 
cntropic forces and consistent with previous experimental findings that nanometer-scale layer widths 
may only be observed very close to the melting point. 
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The term premelting refers to the formation of a thin, 
thermo dynamically stable, liquid-like film at an interface 
for temperatures below the equilibrium melting point 
(Tj\/). Premelting at grain boundaries (GB) can have dra- 
matic consequences in the context of materials process- 
ing, and the physical properties of polycrystals at high 
homologous temperatures. Despite the importance of this 
phenomenon, direct experimental observations of GB 
premelting remain relatively rare, [3,0, S| particularly in 
the case of pure materials 0, Q ■ Consequently, outstand- 
ing fundamental questions remain concerning the nature 
of the forces which drive premelting at these internal in- 
terfaces. In this Letter we introduce a molecular dynam- 
ics (MD) method that exploits large fluctuations in GB 
width to compute short-range forces resulting from the 
overlap of spatially diffuse crystal-melt interfaces from 
two grains of different orientations. We demonstrate the 
application of this method in a direct calculation of the 
excess free-energy of the GB as a function of this width 
for a high-energy boundary in a classical model of el- 
emental Ni. The results yield quantitative insights into 
the relative magnitudes of these short-range structural 
forces and other long-ranged contributions, and help ex- 
plain the origins of the experimental observations 0] that 
GB premelting in pure metals may occur only over ex- 
tremely small temperature ranges near Tm- 

Premelting generally reflects a competition between 
opposing bulk and interfacial thermodynamic factors, 
giving rise to a free energy (per unit area) of the fol- 
lowing form (e.g., [J]): 



G(w) = AGfw + *(«;). 



(1) 



In Eq. [T] w represents the width of the premelted layer, 
AG / is the free energy difference between solid and liquid 
(per unit volume) that penalizes the formation of liquid 
films below Tm, and ^>(w) is the so-called "disjoining po- 
tential" which takes the limits of 7gb (the interfacial free 
energy of a "dry" grain boundary) and (twice the 



solid-liquid interfacial free energy) for zero and infinite 
w, respectively. In general, the disjoining potential con- 
tains both repulsive and attractive contributions. Long- 
ranged dispersion forces lead to an attractive interaction 
between solid-liquid interfaces 0, 0] which are dominant 
at large w and are predicted to give rise to finite inter- 
facial widths at T M 4]. For A7 = j GB — 2"fsL > 0, a 
repulsive contribution to ^(r) arises from short-ranged 
structural interactions (f sr ), associated with the over- 
lap of the diffuse regions of the solid-liquid interfaces. 
The exact nature of this structural contribution remains 
less well understood. 

Mean-field arguments , as well as lattice-gas models 
(e.g., Q), yield an exponentially decaying form for the 
short-ranged contribution to the disjoining potential: 



V sr (w) = 2~/ sl + A 7 exp[-w/<5]) 



(2) 



where S is an interaction length on the order of the 
atomic spacing. In the absence of long-ranged disper- 
sion forces, and neglecting capillary fluctuations which 
lead to an additional repulsive contribution to (see 
below) , insertion of Eq. [2] into Eq. [1] leads to the predic- 
tion of a continuous premelting transition with an equi- 
librium grain boundary width that diverges logarithmi- 
cally as Tm is approached from below. Recent theoreti- 
cal results suggest that the nature of can be much 
more complex. Diffuse-interface theories, [8, 9] which ne- 
glect long-ranged forces and capillary fluctuations, have 
shown that the dependence of w on temperature may in 
some cases display a discontinous jump, with the coex- 
istence of "wet" and "dry" interface states, while other 
parameter choices lead to continuous increases in w up 
to Tm- In recent applications of the phase-field crystal 
(PFC) method to the study of grain boundary premelt- 
ing [13, EH , results for 2D hexagonal systems give a dis- 
joining potential that is purely repulsive above a critical 
misorientation [Iol | , but exhibits a minimum, correspond- 
ing to a finite layer width at the melting point, below it. 
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While these theoretical models thus suggest a rich be- 
havior for the disjoining potential in general, it remains 
unclear for which GB misorientations the various qualita- 
tive forms for ty sr (w) may be expected in real materials, 
and how structural short-range forces compare qualita- 
tively to long range forces. To facilitate further progress 
in the understanding of the forces that drive premclting, 
we describe in the remainder of this Letter a quantitative 
framework for the direct calculation of s S sr (w), through 
histogram analyses of interface widths derived from MD. 

Classical MD simulations provide a framework ideally 
suited for probing the short-ranged structural contribu- 
tions to ^(w). Such simulations have been employed ex- 
tensively in the past to study GB premelting 3, hj and 
in the present work we propose a methodology to extend 
the analysis of such MD results as a framework for ex- 
tracting 9(w), We demonstrate the approach for a clas- 
sical model of elemental Ni, described by the embedded- 
atom potential of Foiles, Baskes and Daw [l6j]. The po- 
tential was chosen as we have previously calculated the 
solid-liquid intcrfacial free energies, melting temperature 
and solid- liquid thermodynamic properties with high pre- 
cision. A value for jsl of 285 mJ/m 2 for the potential has 
been determined using the capillary fluctuation method 
flij ]. and a coexistence technique was used to compute a 
melting temperature of 1710 ± 5K [ijj (from subsequent 
coexistence runs and an analysis of GB width fluctua- 
tion data presented below the uncertainty in this esti- 
mate has been reduced to approximately 1 degree). We 
began by considering a total of four boundaries with a 
range of zero-temperature grain-boundary energies span- 
ning 450 to il43i mJ/m 2 , which is 0.9 to 2.5 times the 
value of 7sl- For each GB we performed a conjugate 
gradient minimization (exploring also the microscopic 
translational degrees of freedom and the excess number 
of atoms at the grain boundary) to derive an optimized 
zero-temperature interface structure. With this structure 
as a starting point the GBs were heated gradually up to 
the melting point employing constant-temperature MD 
simulations following the procedure described in the on- 
line supplemental material. 

Figure [T] shows the calculated excess volume of each 
GB, displaying three qualitatively different behaviors. 
The highest energy boundaries feature an excess volume 
displaying a logarithmic divergence characteristic of a 
continuous premelting transition. The lower energy GBs 
show two different behaviors; in one case the excess vol- 
ume rises with increaing temperature but then plateaus, 
maintaining a finite excess volume at the melting tem- 
perature. The lowest energy GB shows an excess volume 
that is relatively small and only weakly dependent on 
temperature. The range of behavior demonstrated by the 
different GBs in Fig.[T]is qualitatively similar to that seen 
in very recent GB simulations for Si [3] • 

For the remainder of this Letter we focus on one of the 
two high-energy boundaries displaying clear premelting 
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FIG. 1: Calculated excess volumes versus temperature for 
four grain boundaries in elemental Ni. 



behavior, namely the £9 boundary characterized by a 
120° rotation about the GB normal lying along the [511] 
crystallographic direction. We study the structural prop- 
erties of the grain boundary at five temperatures over a 
range of undercoolings 30 to 2 K below Tm- After 4.2 ns 
equlibration runs at each temperature, statistics were ob- 
tained for ro at a given temperature as follows. For each 
snapshot (selected at a frequency of 10 ps) w is deter- 
mined by utilizing the scheme developed in the capillary 
fluctuation method [15| . Each atom is assigned a struc- 
tural order parameter, </>i, constructed from the positions 
of the 12 nearest neighbor atoms and the 4>i values are 
then averaged in bins along the direction normal to the 
boundary. The point of inflection in the average order 
parameter profile is taken as the position of one of the 
solid-liquid interfaces. As described in more detail in the 
accompanying on-line material, the procedure is repeated 
to locate the second solid-liquid interface and hence the 
GB width. 

An important observation in the present work is that 
the width of the GB regions is highly dynamic in the 
MD simulations, particularly at the temperatures closest 
to Tm- This point is illustrated clearly in Fig. [2] which 
shows three snapshots, taken from a 40 ns simulation 
at an undercooling of 2 K, where the atoms have been 
color coded based on their 4n values; blue representing 
a liquid-like environment and red the crystal. The snap- 
shots clearly demonstrate the presence of large fluctua- 
tions in the width of the premeltcd layer over the course 
of the simulation. The highly dynamic nature of the pre- 
melted layer provides a framework for extracting the dis- 
joining potential. We show in Fig. [3] histograms of inter- 
face width obtained at five temperatures near the melting 



FIG. 2: Snapshots from an MD simulation at an undercooling 
of 2 K illustrating the dynamic nature of the GB width. The 
red (blue) areas indicate regions of solid (liquid) like liquid or- 
der. The liquid- like retions at the far left and right hand sides 
represent premelting of the free surfaces of the simulation cell, 
while that in the middle corresponds to the premelted grain 
boundary. 

point. The solid lines represent least-squares fits to the 
data employing the thermodynamic model of Eq. [1] as 
follows. The probability (P(w)) of observing a premelted 
layer width w is given as: 

P(w) = Cexp[~AG(w)/k B T] (3) 

where C is a temperature-dependent normalization con- 
stant, A is the cross-sectional area and G(w) is defined 
in Eq.Q] The data in Fig. Q] suggests a logarithmic diver- 
gence of w with increasing temperature, and in order to 
fit Eq.[3]to the MD data, we therefore employ the form for 
the disjoining potential given in Eq. [5] The least square 
fits of P(w) versus w for all five undercoolings studied 
are shown in Fig. [3] The excellent agreement suggests 
that the free energy of Eq. [1] together with the disjoin- 
ing potential (Eq. [5J represents an accurate model for 
the premelting behavior of this boundary. The analysis 
of Fig. [3] assumed a melting point of Tm = 1710-fT. If in- 
stead a value of just one degree different, i.e., Tm = 1709, 
is assumed, then a poor fit to P(w) is obtained at the 
lowest undercooling. In addition, for simulations run at 
a temperature of YJ12K the system exhibited a gradual 
melting. These findings, together with the results of sep- 
arate coexistence simulations, indicate that the melting 
temperature is known to a precision approaching ±1°. 

As an independent check on the validity of Eq. [2l we 
employ an additional analysis of the data of Fig. [3l From 
Eq. [3] the disjoining potential can be written in terms of 
P{w) as: 

*(«>) = -{k B T/A)lnP{w,Ti) - AG f w + a t (4) 

where the at are unknown constants related to C in Eq. [3] 
and the subscript i denotes a separate histogram of data 
corresponding to each undercooling. The a, can be de- 
termined by a least square fitting procedure such that all 
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FIG. 3: The distribution function P(w) vs w from the MD 
simulations (symbols) versus the least square fits of the pre- 
melting model of Eqs. [T] and [2] 

the data sets can be merged and the entire function *B(w) 
constructed. Notice the procedure adopted here is analo- 
gous to the histogram method, often employed in Monte 
Carlo simulations to extract transition states and energy 
barriers, but with the undercooling playing the role of a 
bias potential. The results of the histogram procedure 
are shown in Fig. [3] The inset of the figure plots the 
right-hand side of Eq. |4j with all the constant terms o» 
set to zero to illustrate that different undercoolings sam- 
ple a range of w regimes of ^(w). The main figure shows 
the final fy(w) function along with a fit (solid line) to 
the exponential form given in Eq. [2] It is important to 
note that the fit parameters obtained via the histogram 
method (S = 2A9A and A7 = 156m J/m 2 ) compare very 
well to those derived through the individual fits of the 
separate histograms in Fig. [3} S = 2.67 ± 0.18A and 
A7 = 127 ± 26mJ/m 2 . 

The present simulations have quantified only the short- 
ranged contribution to ty(w). However, for setting the 
temperature scale over which premelting can be observed 
in elemental metals, we estimate that this is the dominant 
contribution. For values of w ~ 1 nm we can estimate 
an upper bound on the dispersion forces using a value 
of the Hamaker constant measured for surface premelt- 
ing of metals [lj|. This gives a contribution to \&(w=l 
nm) ~ .4 mJ/m 2 , approximately one order of magnitude 
smaller than the results presented in Fig. 4, computed 
from the value of ^> sr using the parameters derived above. 
There is also a known entropic contribution to ^(w) as- 
sociated with the long- wavelen gth fluctuations of the two 
interfaces for large separation [19(. This contribution is 
generally repulsive because the conditions that the two 
interfaces do not intersect reduces the available config- 
urational entropy of interface meandering. While this 
entropy reduction produces a physically important long- 
range force for one-dimensional interfaces, such as step 
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FIG. 4: An illustration of the histogram method used to ex- 
tract the disjoining potential. The inset plots the right hand 
side of Eq.[4]with the constants a,i set to zero. The main plot 
merges the individual histograms of data and reproduces the 
complete disjoining potential 9(w). The solid line is a best fit 
to the exponential decay given in Eq. [2] 



ledges on surfaces, it produces only a subdominant short- 
range forces for two-dimensional interfaces owing to the 
slow logarithmic growth of the mean-square fluctuation 
amplitude with interface area, as compared to the much 
faster square-root growth of this amplitude with interface 
length in one dimension. In particular, a straightforward 
estimate of this force using analytical results of the liter- 
ature [3| and the parameters for Ni from the MD results 
show that this entropic force is completely negligible in 
comparison to the short-range structural forces computed 
here. 

Thus, for high-energy boundaries it can be expected 
that the lowest temperatures where premelting will be- 
come appreciable is set through the relation (AT/Tm) = 
(A'y/LpS) exp (—w eq /5), with w eq ~ 1 nm, where we have 
expressed AG/ = LAT p/Tm in terms of the latent heat 
per atom (L), the solid density (p) and the undercooling 
(AT = Tm—T). The present results give 8 on the order of 
an interatomic spacing and A7 on the order of half 7sl- 
While the exact values will vary somewhat depending on 
the system, we believe these values are well representa- 
tive of high-energy boundaries in pure metals. With these 
estimates we obtain a value of AT required to obtain 
w ~ lnm of (AT/Tm) ~ (ck/2) exp (—4), where we have 
used pS 1 / 3 ~ 1 and the relation j s lp~ 2 / 3 / L — a, where a 
is the Turnbull coefficient [Hj which has a roughly con- 
stant value of about 0.5 for elemental metals [15j. The 
estimate of the undercooling required for a 1 nm pre- 
melted film is thus AT/T M ~ 0.005. The results are 
consistent with the experimental studies of Balluffi and 
co-workers [H who estimated a lower bound of T=0.999 
Tm for the temperature where boundary widths of a few 
nm could be observed experimentally. 
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